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Abstract. Multiscale dynamics are ubiquitous in applications of modern science. Because 

of time scale separation between relatively small set of slowly evolving variables and (typically) 

much larger set of rapidly changing variables, direct numerical simulations of such systems often 

require relatively small time discretization step to resolve fast dynamics, which, in turn, increases 

computational expense. As a result, it became a popular approach in applications to develop a closed 

approximate model for slow variables alone, which both effectively reduces the dimension of the phase 

5_^ ' space of dynamics, as well as allows for a longer time discretization step. In this work we develop a new 

O | method for approximate reduced model, based on the linear fluctuation-dissipation theorem applied 

^•^ , to statistical states of the fast variables. The method is suitable for situations with quadratically 

nonlinear and multiplicative coupling. We show that, with complex quadratically nonlinear and 

Q\ ' multiplicative coupling in both slow and fast variables, this method produces comparable statistics 

to what is exhibited by an original multiscale model. In contrast, it is observed that the results 

from the simplified closed model with a constant coupling term parameterization are consistently 
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1. Introduction. Multiscale dynamics are ubiquitous in applications of mod- 
ern science, with geophysical science and climate change prediction being well-known 
examples [TU1HHEDCI3] • Direct numerical simulations of multiscale systems are dif- 
ficult, both because a relatively short time discretization step is required to resolve 

w~] ■ fast dynamics, and due to a large number of dynamical variables. Moreover, in some 

applications such as climate change prediction one is interested in long-term statistics 
of slow dynamics, which further increases computational expense. 

rf\ ' A popular approach for simulating multiscale dynamics in practice with limited 

computational resources is to create an approximate reduced model for slow variables 
alone, which allows to increase the length of the time discretization step and reduced 
the dimension of the phase space of the system, which is accomplished via an approx- 
imate closure of the coupling terms between slow and fast variables of the system. 
Many closure methods were designed for multiscale dynamical systems [Tll[T3"I[lT?l - r2"2"] , 
based on the averaging formalism for the fast variables 24,29.30 . Some methods 
approximate the coupling terms with appropriate stochastic processes [19 22,131] or 
conditional Markov chains [11], while others [13] parameterize slow- fast interactions 
by direct tabulation and curve fitting. 

In a recent work [3] the author developed a simple approach of computing the 
reduced model for slow variables alone via a single computation of relevant statistics 
for the fast dynamics with a fixed reference state of the slow variables, using the 
linear fluctuation-dissipation theorem (FDT) [1-3.5 9.18,25]. It was shown that, for 
the appropriately rescaled two-scale Lorenz 96 model [5] with linear coupling between 
slow and fast variables, the method reproduced statistics of the slow variables of the 
complete two-scale Lorenz model with good precision. However, nonlinear coupling 
was not addressed in pQ. 
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This work is a natural extension of the method in [1] onto quadratically nonlinear 
and multiplicative types of coupling in both slow and fast variables. It is also based 
on the linear FDT, which is used to compute the response of linear, nonlinear and 
multiplicative terms of the fast variables coupling to changes in the slow variables. 
We show numerical experiments for this method with the two-scale Lorenz model and 
general forms of coupling which include multiplicative and quadratically nonlinear 
terms. 

The manuscript is organized as follows. In Section [5] we formulate the theory 
for the new method. In Section [3] we present the two-scale Lorenz model [T51[TB1[T7] . 
rescaled in such a way that the mean states and variances of both the fast and slow 
variables are near zero and one, respectively [4j[5]. Section [4] shows the results of 
numerical simulations with both the two-scale Lorenz model and the reduced model 
for slow variables only, comparing different statistics of the time series. Section [5] 
summarizes the results of this work. 

2. Derivation of the reduced model. We start with a general two-scale sys- 
tem of differential equations of the form 

^=F(x,y), ^ = G(x,y), (2.1) 

where x — x(t) G TBL Nx are the slow variables, y = y(t) 6 R^" are the fast variables, 
and F and G are N x and N y vector- valued functions of x and y, respectively. We 
assume that the y-variables are sufficiently fast for a valid approximation of the system 
in (|2.ip by the averaged dynamics for a;, given by 

ri-r f 

= (F)(x) } (F)(x)= F(x,z)d^(z), (2.2) 

M J m N y 

for finite times (for a more detailed description of the averaging formalism, see [2J[51 
[MHHllSn])- Here, /i x denotes the invariant probability measure of the limiting fast 
dynamics, which are given by system 

^ = G(x,z). (2.3) 



Above, s is a constant parameter, and the solution of (|2.3p is given by the flow 
z(t) — 4> x zq. We tacitly assume that all typical initial conditions Zq fall into the 
support of the same ergodic component of \i Xl and that (F)(x) varies smoothly with 
respect to x, as it often happens when n x is an SRB measure [12,26 28,32 . Using 
the ergodicity assumption for fj, x , we can practically compute the measure average 
via the time average 



1 r 

(F)(x)= lim - / F(x,z(r))dT, 
r->-oo r J 



(2.4) 
'" Jo 

where z(t) is a long-term trajectory of (I2.3[) . 

Following [4], here we propose an approximation to (F)(x), based on the linear 
FDT. In order to do this, certain assumptions have to be made regarding coupling, 
that is the dependence of F and G on the fast and slow variables, respectively. In [4] 
it was assumed that F depends linearly on the fast variables, and G depends linearly 
on the slow variables (linear coupling). Here we assume that F is quadratic in the 
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fast variables in the vicinity of z(x), which is the mean state of (|2.3|) with x set as a 
constant parameter: 

OF 
F(x, z) = F(x, z{x)) + j-ix, z{x)){z - z{x))+ 

uf- F V M 

where dF/dy and d 2 F/dy 2 denote the first and second partial derivatives of F with 
respect to its second argument, respectively, and ":" is the element-wise (Hadamard) 
matrix product with summation. The assumption of quadratic dependence of F 
on the fast variables is not too restrictive for practical applications; indeed, many 
real-world geophysical processes have at most quadratic dependence on velocity or 
streamfunction fields due to advection. Now, the average with respect to [i x is given 

by 

ld 2 F . 
2W 



(F)(x) = F(x, z(x)) + - — (*, z{x)) : £(*), (2.6) 



where £(2;) is the covariance of z, centered at z(x). The mean state z{x) and the 
covariance S(a;) are given by 

z(x) = I zdfi x (z) 7 (2.7a) 

Jr n v 

£(cc) = / {z-z{x))®(z-z(x))dn x (z). (2.7b) 

Jr n v 

As we can see, the average of F(x, z) with respect to /i x is now expressed as a 
nonlinear function in the mean state z(x), and linear function in covariance H(x). 
If we know how these quantities respond to changes in x, we can also calculate the 
approximation of the ^-dependent average (F)(x). 

Now, the response of z(x) and E(sc) to changes in x can be estimated via the 
linear FDT. In order to obtain the response formulas, we have to impose some restric- 
tions on the structure of G(x, y) in y. Here, we assume that G(x, y) can be written 
as 

G(x,y)=g(y) + H(x)y + h(x), (2.8) 

where g(y) is a iV y -vector nonlinear function of y, h(x) is a A^-vector function of 
x, and H(y) is a N y x TVy-matrix valued function of x. Again, this assumption is 
not too restrictive for practical applications, as only the nonlinear part of G does not 
depend on x. In this case, the limiting system in (J2.3I) can be written as 

P=g(z) + H(x)z + h(x), (2.9) 

where a; is a constant parameter. Now, let H* = (H(x)), h* = (h(x)) denote the 
long-term averages of H(x) and h(x) over a trajectory of (|2.1|) . so that we can write 
the fast limiting system as 

— = g(z) + (H* + 5H{x))z + {h* + 6h(x)), 
SH(x)=H(x)-H*, Sh(x) = h(x) - h* . 
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Given the mean state z* and the mean-centered covariance matrix S* for unperturbed 
(|2.10p with SH(x) and 8h(x) set to zeros, one can think of z(x) and S(a;) as responses 
of z* and S* to nonzero 5H(x) and 8h(x). These responses can be written as linear 
approximations via the FDT: 



z(x) w z* + R h ^ z (Sh(x) + 8H{x)z*) + R H ^ z SH(x), 



H^fT,i 



E(sc) « S* + iT^tfhix) + SH(x)z*) + R"^5H{x), 



(2.11) 



where the linear response operators R l z , R z , R , R are given by the 



quasi-Gaussian formulas 



R 



h^>z 



lim - / (Zi(t + s)- z*){z k {t) - z* k ) di 



dszij 1 , 



(2.12a) 



T/H^-z _ 
n ijk — 


/•OO 

Jo 


i r 

lim - / (z;(i + s) - z*)(zi( 

r^oo r J 


t)-z?)> 




x(z k (t) -z* k )&t 


ds S- 1 


/>oo 


1 r 

lim - / (zi(t + s) - zf)(zj(t + s) - z*)x 

r—¥oo v Jo 


n ijk — 

Jo 




x{zi(t)-zf)dt 


ds s r fe s 


/•OO 

^ijkl = / 

Jo 


lim 

r— >oo 


i r 

- / (z J (£ + s)-z*)(z J (i + 
r Jo 


s 


)-f*)x 



x(z m (i) - O (**(*) - **) di Knk - S *A 



ds. 



(2.12b) 



(2.12c) 



(2.12d) 



Above, the time averaging is performed over a single long-term trajectory of the 
unperturbed fast limiting system in (|2.10p . with 8h and 6H set to zero. The details of 
derivation of (|2. 1 1|) and (j2.12[) are given in Appendix^ along with relevant references. 
Combining (|2.11[) with (|2.6[) yields the approximate reduced system for slow vari- 
ables alone. Indeed, observe that if the averages z* , £*, and the response operators 
R h ^ z , R H ^ Z , R h ^, and R H ^ are computed, then (12TTTJ) and, therefore, flU) are 
known explicitly for given parameter x. 

3. The Lorenz model with nonlinear coupling. The rescaled Lorenz model 
with nonlinear coupling is given by 



1 F x — 

Xi-l(Xi+l - Xi- 2 ) + —{x(x i+ i - Xi-2) - X%) + X o9 



A 



J2 [(a + b Xi )y hJ + (c + dXiXylj - 1)] 



(3.1a) 



ViJ 



1 p _ y 

Vi,j+i(Vi,j-i - Ui.3+2) + -z-(v(Vi,j-i ~ yij+2) ~ Vi,j) + 



I3y 
+ —[(« + cy%,j)xi + (b + dyij)(xf - 1)] , 



ft 



(3.1b) 
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with 1 < i < N x and 1 < j < J, such that N y = N X J . The model has periodic 

boundary conditions: x i+ ^ x — Xi, Vij+j = Vi+ij, and 2/i+A^j = Vi,j- The parameters 
(x,[3 x ) and (y,(3y) are the (mean, standard deviation) pairs for the corresponding 
uncoupled and unrescaled Lorenz models 

Xi = Xi-i(x i+ i - Xi- 2 ) - Xi + F x , (3.2a) 

Vi,j = yi, j+ i{yi,j-i - yt. 3 +2) - Vij + F y , (3.2b) 

with the same periodic boundary conditions. The rescaling above ensures that the 
Lorenz model in (j3.ll) has zero mean state and unit standard deviation for both slow 
and fast variables in the absence of coupling (X x = X y = 0), and remain near these 
values when A^ and X y are nonzero. The original rescaled Lorenz model in [4,5 with 
linear coupling corresponds to the set of parameters a = 1, b = c = d = 0. The 
nonlinear coupling above preserves the energy of the form 

, N x . N x J 

* = ^£*? + ^££&- (3-3) 

Indeed, observe that 

N NT 

&E \ x £^ . X 



/ j x i x i + j 2.^1 2.^1 yi^yi'i 



2—1 2—1 J — l 

A A N °> J 
= - -^j- £ £ (aXiVij + bxfy itj + cx.yfj + dx 2 y? tj ) + (3.4) 

i=l j=l 

A A N " J 
+ ~^j- £ £ (aXiyij + bx 2 y tJ + cx l y 2 ij + dxjy 2 ^) = 0. 

t=l j=l 



At this point, one can see that h(x) and H(x) are given by 

hij(x) = — (aXi + bxf) , 



e 



Hi,j,i',j'(y) = -f5\f (cx l +dx 2 ) 



(3.5) 



in particular, H(x) is a diagonal matrix. Also, 

i0(*) : V(x) = -^ (c + d Xi )J2^(x), (3.6) 

that is, only the diagonal entries of the covariance matrix are needed. This results in 
R h ^ s , R h ^, R H ^ S and R H ^ all being matrices rather than 3- or 4-dimensional 
tensors. 

4. Numerical simulations. Here we show the results of numerical simulations 
with the new reduced model for slow dynamics of the rescaled Lorenz model with 
nonlinear coupling in (|3.1[) . In particular compare the numerical simulations for the 
three following systems: 
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1. The full two-scale rescaled Lorenz system from (j3.ll) ; 

2. The reduced model for slow dynamics from (|2.6|) and (|2.11l) ; 

3. A simplified version of (12.61) with both the mean state z and covariance matrix 
£ of the fast variables are fixed at their average values z* and £*, without the 
correction terms from (|2.11[) (further referred to as the "zero-order" system). 

The quasi-Gaussian approximations in (|2.12j) are used to compute the approximations 
of the coupling terms in (|2.11|) . For the time averaging in (|2.12j) we use long-term 
trajectories with averaging time window equals 10000 time units, while the correlation 
time window equals 50 time units (it was observed that the time autocorrelation 
functions in (|2.12[) decays essentially to zero within the 50 time-unit window for all 
studied regimes). The reference mean state z* and the covariance matrix S* are 
also computed by time-averaging of the full two-scale Lorenz model with the same 
averaging window of 10000 time units. 

The statistics of the model are invariant with respect to the index permutation 
for the variables Xi, due to translational invariance of the studied models. For the 
numerical study, we compute the following long-term statistical quantities of Xt: 

a. The probability density functions (PDF), computed by standard bin-counting. 

b. The time autocorrelation functions (xi(t)xi(t + s)), where the angled brackets 
denote the time average is over t. These autocorrelation functions are normalized 
by the variance (xf), so that the initial value at s — is always 1. 

c. The time cross-correlation functions (xi(t)Xi-\-\(t + s)), also normalized by the 
variance (xf). 

d. The energy autocorrelation function 

{xUt)x1{t + s)) 



K{s) = 



(x 2 i ) 2 +2(x i (t)x i (t + s))^ 

This energy autocorrelation function measures the non-Gaussianity of the process. 

It is identically 1 for all s if the process is Gaussian (the Ornstein-Uhlenbeck process 

being an example). For details, see |21j . 
The performance of the proposed approximation of the slow dynamics depends on 
several factors. First, the precision will be affected by the non-Gaussianity of the 
fast dynamics, since the quasi-Gaussian linear response approximation is used for the 
computation of the response operators. Second, the dependence of the mean state 
z(x) and covariance E(x) for the fast variables depends on the perturbations 8h(x) 
and SH(x) is generally nonlinear, and that should also affect the precision of the 
approximation. Here we study the behavior of the proposed approximation in variety 
of dynamical regimes of the rescaled Lorenz model in (13.11) . The following dynamical 
regimes are studied: 

• N x = 20, J = 4 (so that N y = 80), so that the number of the fast variables 
is four times greater than the number of the slow variables. 

• £ = 0.01. Typical geophysical processes, such as the annual and diurnal 
cycles, have the time scale separation of roughly two orders of magnitude. 

• F x — 6, F y — 12. The slow forcing F x adjusts the chaos and mixing properties 
of the slow variables, and in this work it is set to a weakly chaotic regime 
F x = 6. The fast forcing F y regulates chaos and mixing at the fast variables, 
which are usually more chaotic and mixing than the slow variables, so it is 
set to F y = 12. 

• A x = Xy = 0.3. This value of the coupling constant is chosen based on the 
previous work [3j , where the same value was used to test the method for linear 
coupling. 
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Fig. 4.1. Probability density functions, time autocorrelations, cross-correlations and energy 
autocorrelations. Parameters: (a, b, c, d) = (0, 0.8, 0, 0). 



• We test a mixture of different constants {a,b,c,d) for each coupling term. 
First, we test the isolated coupling regimes with all constants but one set to 
zero, and then show the results for mixed coupled regimes, where different 
types of coupling are present simultaneously. 

4.1. Single type coupling. In this section we test the regimes with single type 
coupling (one of the constants a,b,c,d is nonzero, and the rest are zero). In this 
section, we only test the regimes with nonzero constants 6, c, d, as the single type 
coupling regime with a ^ corresponds to the model previously studied in [4] . 

4.1.1. Regime with (a,b,c,d) — (0, ±0.8, 0, 0). Here we test the regimes with 
b = ±0.8, and a = c = d = 0. This regime corresponds to bilinear multiplicative 
coupling in the slow variables, and quadratic additive coupling in the fast variables. 
The probability density functions, time autocorrelation, cross-correlation, and energy 
correlation functions are shown in Figures |4 . 1 1 and |4 . 2 1 while the L2-errors between the 
curves are shown in Table |4~T] There is little difference between the plots, which means 
that the fast mean state z(x) and S(x) are not very sensitive to changes in x for 
this type of coupling. Yet, one can see that the reduced model with linear correction 
for z(x) and S(x) more precisely captures statistics of the full scale model, than the 
zero-order reduced model with fast mean state and covariance fixed at z* and S* (see 
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Fig. 4.2. Probability density functions, time autocorrelations, cross-correlations and energy 
autocorrelations. Parameters: (a, 6,c, d) = (0, — 0.8, 0, 0). 
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Table 4.1 

L2-errors between the statistics of the slow variables of the full two-scale Lorenz model and the 
two reduced models, with (a,b,c,d) = (0, ±0.8, 0, 0). Notations: "Reduced" stands for the reduced 
model from 112.61 1 and 1 12. lift , and "Zero-order" stands for the poor man's version of the reduced 
model, with linear approximations for z(x) and S(a;) replaced by constant mean values z* and S*. 



Table FTTj) . Additionally, it appears that there is a weak effect of chaos suppression 
for this type of coupling regardless of the sign of b (for both negative and positive b, 
the statistics of the coupled model show somewhat slower decay of autocorrelations 
and cross-correlations, and more sub-Gaussian energy correlation functions). 

4.1.2. Regime with (a,b,c,d) = (0, 0, ±0.3,0). Here we test the regimes with 
c = ±0.3, and a = b = d = 0. This regime corresponds to quadratic additive 



A closure for slow dynamics of a multiscale system: nonlinear coupling 



N x =20, N =80, ^=0.3, ^=0.3, F x =6, F =12 
(a,b,c,d)=(0,0,0.3,0), 8=0.01 



N x =20, N =80, ?i x =0.3, A. =0.3, F x =6, F =12 
(a,b,c,d)=(0,0,0.3,0), e=0.01 




N x =20, N =80, A. x =0.3, X =0.3, F x =6, F =12 
(a,b,c,d)=(0,0,0.3,0), e=0.01 



Time 

l.3,X 

(a',b,c,d)=(0,0,0.3,0), e=0.01 



N„=20, N y =80, ?i y =0.3 



^=0.3, F x =6, F - 

.3,r "-' ' 



12 



Full system 

Reduced system 

Zero-order system 




Full system 

Reduced system 

Zero-order system 




Time 



Time 



Fig. 4.3. Probability density functions, time autocorrelations, cross-correlations and energy 
autocorrelations. Parameters: (a, 6, c, d) = (0,0, 0.3, 0). 



coupling in the slow variables, and bilinear multiplicative coupling in the fast variables. 
The probability density functions, time autocorrelation, cross-correlation, and energy 
correlation functions are shown in Figures I4.3I and 14.41 while the L2-errors between 
the curves are shown in Table |4~21 Like in the previous case, there is little difference 
between the plots, which means that the fast mean state z(x) and £(a;) are not very 
sensitive to changes in x for this type of coupling. Yet, one can see that the reduced 
model with linear correction for z(x) and £(a:) more precisely captures statistics of 
the full scale model, than the zero-order reduced model with fast mean state and 
covariance fixed at z* and X* (see Table |4T^1) . Additionally, it appears that there is a 
weak effect of chaos suppression for this type of coupling for positive value of c = 0.3 
(the statistics of the coupled model show somewhat slower decay of autocorrelations 
and cross-correlations, and more sub-Gaussian energy correlation). No such effect is 
observed for the negative value c = —0.3. 

4.1.3. Regime with (a,b,c,d) — (0, 0, 0, ±0.3). Here we test the regimes with 
c = ±0.3, and a = b = d = 0. This regime corresponds to quadratic multiplica- 
tive coupling in both the slow and fast variables. The probability density functions, 
time autocorrelation, cross-correlation, and energy correlation functions are shown in 
Figures 14.51 and 14.61 while the L2-errors between the curves are shown in Table 14.31 
Unlike previous cases, there is a significant difference between the results of the full 
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Fig. 4.4. Probability density functions, time autocorrelations, cross-correlations and energy 
autocorrelations. Parameters: (a, b, c, d) = (0,0, — 0.3, 0). 
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= 0.3 




(a 


b,c,d) = (0,0,0.3,0) 


(a, 


b,c,d) = (0,0, 


-0.3,0) 




Reduced 


Zero-order 




Reduced 


Zero-order 


PDF 


5.49 -10- 3 


2.314 -10- 2 


PDF 


3.707 -10- 3 


1.447 • 10~ 2 


Corr. 


5.559 -10- 2 


0.1788 


Corr. 


5.238 -10- 2 


0.2399 


C-corr. 


6.265- 10~ 2 


0.2055 


C-corr. 


6.271- 10~ 2 


0.2845 


K-corr. 


9.605 • 10" 3 


2.659 -10" 2 


K-corr. 


9.62 -10~ 3 


2.041 • 10~ 2 



Table 4.2 

L,2-errors between the statistics of the slow variables of the full two-scale Lorenz model and the 

two reduced models, with (a,b,c,d) = (0, 0, ±0.3, 0). Notations: "Reduced" stands for the reduced 

model from 112.61 1 and 1 12. lit , and "Zero-order" stands for the poor man's version of the reduced 

model, with linear approximations for z(x) and S(a;) replaced by constant mean values z* and X*. 



two-scale/first-order reduced model, and the zero-order reduced model, which means 
that the fast mean state z(x) and S(a;) are sensitive to changes in x for this type 
of coupling. Similar to the previous cases, one can see that the reduced model with 
linear correction for z{x) and Tl(x) more precisely captures the statistics of the full 
two-scale model, than the zero-order reduced model with fast mean state and co- 
variance fixed at z* and X* (see Table |4~3"| . Unlike previous types of coupling, here 
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N x =20, N =80, A, x =0.3, A. =0.3, F x =6, F =12 
(a,b,c,d)=(0,0,0,0.3), 8=0.01 



N x =20, N =80, ?i x =0.3, X y =0.3, F x =6, F =12 
(a,b,c,d)=(0,0,0,0.3), e=0.01 




N x =20, N y =80, A.. 



=0.3, X v =0.3, F x =6, F v =12 
',0.3' - """ ' 
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Time 



N x =20, N =80, V - 3 . K=0-3, F x =6, F =12 
(a,b,c,d)=(0,0,0,0.3), e=0.01 



Full system 

Reduced system 

Zero-order system 




Full system 
Reduced system 
Zero-order system 




Time 



Time 



Fig. 4.5. Probability density functions, time autocorrelations, cross-correlations and energy 
autocorrelations. Parameters: (a, b, c, d) = (0, 0, 0, 0.3). 



the chaos suppression or amplification depends on the sign of the constant d. For 
the positive value d = 0.3, the coupled dynamics and the reduced model are clearly 
less chaotic and mixing than the zero-order reduced model with fast mean state and 
covariance fixed at z* and £*, which follows from the difference in decay of the cor- 
relation functions. The opposite effect is observed for the negative value d = —0.3. 
In particular, observe that the PDF of the zero-order reduced model displays three 
peaks for d = —0.3 (which is a sign of quasi-periodic motion), while the PDFs of both 
the two-scale Lorenz model and reduced models are unimodal. 

4.2. Combined coupling. In this section we test the regimes with all types of 
coupling observed previously in Section [4.1[ combined together in the Lorenz model. 
In addition, the linear coupling is also switched on by setting a = 1. Here we study 
the sets of parameters (a,b,c,d) = (1, ±0.8, ±0.3, ±0.3). The probability density 
functions, time autocorrelation, cross-correlation, and energy correlation functions 
are shown in Figures 14.71 and 1 4.81 while the L2-errors between the curves are shown in 
Table [4~4l It turns out that when all types of coupling are combined together there is a 
major difference between the two-scale/first-order reduced models, and the zero-order 
reduced model with fast mean state and covariance fixed at z* and £ , which means 
that the fast mean state z(x) and £(x) are quite sensitive to changes in x. Similar 
to the previous cases, here one can see that the reduced model with linear correction 
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=20, N =80, V - 3 ' V - 3 ' F x= 6 > F v= 12 
(a;b,c,d)=(0,0,0,-0.3), £=0.01 



6, F„=12 



N x =20, N =80, X„=0.3, ?i =0.3, F x -„. , 
(a,b,c,d)=(0,0,0,-0.3),e=0.01 




Time 



Time 



Fig. 4.6. Probability density functions, time autocorrelations, cross-correlations and energy 
autocorrelations. Parameters: (a, b, c, d) = (0, 0, 0, — 0.3). 
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PDF 


1.342- 10~ 3 


5.272 • 10~ 2 


Corr. 


7.397 -10- 2 


0.2638 


Corr. 


1.923- 10" 2 


0.3881 


C-corr. 


9.183 -10- 2 


0.3196 


C-corr. 


2.497- 10~ 2 


0.4624 


K-corr. 


1.632 • 10" 2 


5.002 -lO" 2 


K-corr. 


3.656 • 10" 3 


6.977 -10" 2 



Table 4.3 

L,2-errors between the statistics of the slow variables of the full two-scale Lorenz model and the 

two reduced models, with (a,b,c,d) = (0, 0, 0, ±0.3). Notations: "Reduced" stands for the reduced 

model from 112.61 1 and 1 12. lit , and "Zero-order" stands for the poor man's version of the reduced 

model, with linear approximations for z(x) and S(a;) replaced by constant mean values z* and S*. 



for z(x) and S(x) much more precisely captures statistics of the full two-scale scale 
model, than the zero-order reduced model (see Table l4~3l) . Here the chaos suppression 
or amplification depends on the signs of the constant 6, c and d. For positive values of 
these constants, the coupled dynamics and the reduced model are clearly less chaotic 
and mixing than the zero-order reduced model, which follows from the difference in 
decay of the correlation functions. The opposite effect is observed for negative values 
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N x =20, N =80, X x =0.3, ?i =0.3, F x =6, F =12 
(a,b,c,d)=(1, 0.8,0.3,0.3), £=0.01 



N x =20, N =80, >. x = - 3 . V 03 ' F x= 6 ' F y= 12 
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N x =20, N =80, A x =0.3, X =0.3, F x =6, F =12 
(a,b,c,d)=(1 ,0.8,0.3,0.3), e=0.01 



N x =20, N =80, V - 3 

a,b,c,d)=(1 ,0.8,0.3,0.3), e="0.01 




Time 



Time 



Fig. 4.7. Probability density functions, time autocorrelations, cross-correlations and energy 
autocorrelations. Parameters: (a, b, c, d) = (1, 0.8, 0.3, 0.3). 



of b, c and d. In particular, observe that the PDF of the zero-order reduced model 
displays three strong peaks for d = —0.3 (which is a sign of quasi-periodic motion), 
while the PDFs of both the two-scale Lorenz model and reduced models are unimodal. 
On the contrary, for d = 0.3 the PDF of the two-scale Lorenz model displays three 
peaks, which are roughly captured by the reduced model, while the zero-order reduced 
model produces nearly Gaussian PDF. 

5. Conclusions. In this work we develop a simple approach for approximation 
of multiscale dynamics with nonlinear and multiplicative coupling via a reduced model 
for slow variables alone. The method is based on the linear approximation of averaged 
coupling terms by means of the fluctuation-dissipation theorem, which is only requires 
a single computation of certain long-term statistics from the fast limiting system with 
the slow terms set as constant parameters. This work is a direct extension of [3] 
onto nonlinear and multiplicative coupling (the original method in [3] was developed 
for linear coupling in both slow and fast variables). We verify through the numerical 
simulations with the rescaled two-scale Lorenz 96 model [4][5] that, with nonlinear and 
multiplicative coupling in both slow and fast variables, the new simple reduced model 
produces statistics which are consistent with those of the complete two-scale Lorenz 
model. In contrast, the "zero-order" reduced model with constant parameterization 
of fast variables in coupling terms fails to reproduce the same set of statistics with 
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Fig. 4.8. Probability density functions, time autocorrelations, cross-correlations and energy 
autocorrelations. Parameters: (a, b, c, d) = (1, — 0.8, — 0.3, — 0.3). 
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Corr. 
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C-corr. 
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C-corr. 


0.2477 


1.337 


K-corr. 


0.134 


0.3356 


K-corr. 


2.89- 10- 2 


0.1707 



Table 4.4 

L2-crrors between the statistics of the slow variables of the full two-scale Lorenz model and the 
two reduced models, with (a,b,c,d) = (1, ±0.8, ±0.3, ±0.3). Notations: "Reduced" stands for the 
reduced model from 1 12. 611 and 12.111 . and "Zero-order" stands for the poor man's version of the 
reduced model, with linear approximations for z(x) and S(aj) replaced by constant mean values z* 
and S*. 



comparable precision. The method appears to be convenient for practical applications 
due to its explicit construction - it lacks unknown parameters which have to be 
determined implicitly by comparing the performance of the reduced model against 
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the full multiscale dynamics. 

Appendix A. Detailed derivation of the closure terms for the mean 
state and covariance matrix. 

In order to derive the formulas for the linear responses of the mean and covariance 
of the fast variables, for convenience we first make a linear change of variables in the 
following way: we rewrite (|2.10p as 

^ = S- 1 g(z* + Sq) + S- 1 (H*+5H(x))Sq+ 

dr (A.l) 

+S~ 1 (h* + H*z* + Sh(x) + SH(x)z*), 

where S = X* 5 , and z = z* + Sq, that is, q is the fluctuation of z around z* with 
the identity covariance matrix. Changing notations, we obtain 

-^ = S- l g{z* + Sq) + (H* + 5H{x))q + h + Sh(x), 
dr 

h =S~ 1 (h* +H*z*), H* = S l H*S, ( A - 2 ) 

Sh{x) = S^iShix) + SH(x)z*), SH(x) = S~ 1 dH(x)S. 

Here, we can consider (|A.2[) as a dynamical system perturbed by Sh and SH , which 
has zero mean state and identity covariance matrix in the unperturbed state, and 
use the FDT to estimate the linear response of the mean state (q) and the covariance 
(qq T ) to these perturbations, which can later be mapped into the original coordinates 
via backward linear transformation. The linear responses of the mean state 5{q) and 
covariance S(qq T ) are given by, respectively, 



5(q)i = Rl^Sh, + R?? z 6H jk) 

5{qq T ) l3 = rijjfth + R^SHku 
where the linear response operators are given by 



(A.3) 



d 



R-j I I ^Q-SPQiiMq) d«, (A.4a) 



- * r°° r Ft 

R m Z = / jr^vUkd^ds, (A.4b) 

Jo Jm N y OQj 

R W= / 7 -^ s q®<p s q\ ] d^{q)d S , (A.4c) 

Jo Jwl n v oik 

R "^= / -^^ s q®<P s qh qi d^q)d Sl (A.4d) 

J Jo Jr n v "Ik 

where <fi s is the flow generated by the unperturbed system in (|A.2j) with 5h and SH set 
to zeros, and d/i is its invariant measure [THE]- At this point, we are going to assume 
that jx has a Gaussian density with zero mean state and identity covariance matrix 
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(so-called quasi-Gaussian FDT approximation, 4,6 8,18 ). Then, after integrating 
by parts, the responses of the mean state and covariance can be computed as 



R 



h^rz 



JR"y 



s q) i q j dfi(q) ds, 



(A.5a) 



R 



H^z 
ijk 



s q)iqjqkd/j,(z)ds, 



(A.5b) 



R 



ijk 



/ / {(j) s q®(t) s q) lJ qkd^{z)ds, 

JO JR N v 



(A.5c) 



R 



ijkl 



s q <g> (j> s q)ijqkqi d^{z) - SijS k i 



ds. 



(A.5d) 



By replacing measure averages with time averages (under the ergodicity assumption), 
we obtain 



m 



i 



lim - / qi(t + s)qj(t) dt 



ds, 



(A.6a) 



"■ijk 



lim-/ qi(t + s)qj(t)q k (t)dt 



ds. 



(A.6b) 



R 



ijk 



lim 



q l (t + s)q J (t + s)q k (t)dt 



ds, (A.6c) 



lim - / qi{t + s)q 3 (t + s)q k {t)qi{t) dt - 5ij5 k i 



"'ijkl 



Returning back to the original response coordinates, we obtain 

1 



R> 



lim 

r— >co t 



{z l {t + s)-z*)q J {t)dt 



ds, 



ds. (A.6d) 



(A.7a) 



R 



H^z 

ijk 



lim 



(z l {t + s)-z*)q j {t)q k {t)dt 



ds, 



(A.7b) 



^ijk 



i r 

lim - / (zi (t + s 



)-z*)(zj(t + s)-z*)q k (t)dt 



ds, (A.7c) 



"ijkl 



lim - / ( Zi (t + s) - z*)(zj(t + s) -z*)x 

r— >oo T In 



x« fc (t)gi(t)dt-Ej,.«fcj 



ds. 



(A.7d) 
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Returning back to the original perturbation coordinates, we obtain 



R 



h^z 



lim - I ( Zi (t + s)-zt)(z k (t)-z* k )dt 



d*Kj\ 



(A. 



£Xj ijk 



lim - / (Zi(t + s) - z*)(zi(t) - z*)x 



x(z k {t)-z* k )dt 



<±sZ*r\ 



(A.8b) 



R 



h-yx 
ijk 



lim - / (zi(t + s) - z*)(zj(t + s) — z*,)x 



x( Zl (t)-zf)dt 



te^-\ 



(A.8c) 



R 



ijkl 



lim - / (zi(t + s) - z*)(zj(t + s) - z*)x 



ds, 



x (z m {t) - z*J( Zl (t) - zt) dt E^ 1 - EU k 



(A.8d) 



which is given above in (|2.12[) . 
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